Predicting gene expression from leaf-level spectral reflectance in cotton

Calculate the soft thresholding power (beta)

# Set powers to sample
powers = c(c(1:10), seq(from = 12, to = 60, by = 4))

# Call the network topology analysis function
sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 0)
k <- softConnectivity(datE = datExpr, power = sft$powerEstimate)

# Plot the results:
par(mfrow = c(1, 2))
cex1 = 0.9

# Scale-free topology fit index as a function of the
# soft-thresholding power
plot(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[,
    2], xlab = "Soft Threshold (power)", ylab = "Scale Free Topology Model Fit,signed R^2",
    type = "n", main = paste("Scale independence"))
text(sft$fitIndices[, 1], -sign(sft$fitIndices[, 3]) * sft$fitIndices[,
    2], labels = powers, cex = cex1, col = "red")

# this line corresponds to using an R^2 cut-off of h
abline(h = 0.8, col = "red")

# Mean connectivity as a function of the soft-thresholding
# power
fig <- plot(sft$fitIndices[, 1], sft$fitIndices[, 5], xlab = "Soft Threshold (power)",
    ylab = "Mean Connectivity", type = "n", main = paste("Mean connectivity"))

text(sft$fitIndices[, 1], sft$fitIndices[, 5], labels = powers,
    cex = cex1, col = "red")

par(mfrow = c(1, 2))
hist(k)
scaleFreePlot(k, main = "Check Scale free topology\n")

Blockwise module detection

net = blockwiseModules(datExpr, maxBlockSize = 5000, power = 30,
    TOMType = "unsigned", minModuleSize = 10, reassignThreshold = 0,
    mergeCutHeight = 0.05, numericLabels = TRUE, pamRespectsDendro = FALSE,
    saveTOMs = FALSE, verbose = 0)

# Plot the dendogram with color assignment
mergedColors = labels2colors(net$colors)
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
    "Module colors", dendroLabels = FALSE, hang = 0.03, addGuide = TRUE,
    guideHang = 0.05)

moduleLabels = net$colors
moduleColors = labels2colors(net$colors)
MEs = net$MEs
geneTree = net$dendrograms[[1]]

bwnet = blockwiseModules(datExpr, maxBlockSize = 5000, power = 30,
    TOMType = "unsigned", minModuleSize = 10, reassignThreshold = 0,
    mergeCutHeight = 0.05, numericLabels = TRUE, saveTOMs = FALSE,
    verbose = 0)

# Relabel blockwise modules
bwLabels = matchLabels(bwnet$colors, moduleLabels)

# Labels to colors for plotting
bwModuleColors = labels2colors(bwLabels)

# Comapre the single block to the block-wise network
# analysis
plotDendroAndColors(geneTree, cbind(moduleColors, bwModuleColors),
    c("Single", "Blockwise"), main = "Single block gene dendrogram and module colors",
    dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05)

figure <- plotDendroAndColors(geneTree, bwModuleColors, "", main = "Single block gene dendrogram and module colors",
    dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05)

class(figure)

## [1] "list"

Merging of modules whose expression profiles are very similar

# Merge modules whose expression profiles are very similar
# Calculate Eigengenes to quantify co-expression similarity
# of entire modules
MEList = moduleEigengenes(datExpr, colors = bwModuleColors)
MEs = MEList$eigengenes

# Calculate dissimilarity
MEDiss = 1 - cor(MEs)

# Cluster module eigengnes and Plot thhe results
METree = hclust(as.dist(MEDiss), method = "average")
plot(METree, main = "Clustering of module eigengenes", xlab = "",
    sub = "")

# Create Height Cut
MEDissThres = 0.25
abline(h = MEDissThres, col = "red")

merge = mergeCloseModules(datExpr, bwModuleColors, cutHeight = MEDissThres,
    verbose = 0)
mergedColors = merge$colors
mergedMEs = merge$newMEs

Plot clustereing among samples

# Cluster tress (check if any outliers existed among
# samples)
datExpr_tree <- hclust(dist(datExpr), method = "average")
par(mar = c(0, 5, 2, 0))
plot(datExpr_tree, main = "Sample clustering", sub = "", xlab = "",
    cex.lab = 2, cex.axis = 1, cex.main = 1, cex.lab = 1, cex = 0.5)

# Verify epigenegenes of modules in both approaches are
# extremely similar First calculate the module eigengenes
# for first approches
singleBlockMEs = moduleEigengenes(datExpr, moduleColors)$eigengenes
blockwiseMEs = moduleEigengenes(datExpr, bwModuleColors)$eigengenes

# Now match the eigengens by name and calcualte the
# correlation of the corresponding eigengenes
single2blockwise = match(names(singleBlockMEs), names(blockwiseMEs))

blockwiseMEs

# signif(diag(cor(blockwiseMEs[, single2blockwise],
# singleBlockMEs)), 3)